﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.Statistics
{
    /// <summary>
    /// The model function for fitting a curve.
    /// </summary>
    public delegate GeneralVector ModelFunction(GeneralVector x, GeneralVector parameters);

    /// <summary>
    /// Describes if a parameter should be exclude from the fitting process.
    /// </summary>
    public enum CurveFittingParameter
    {
        /// <summary>
        /// The parameter will be computed by the algorithm.
        /// </summary>
        Free = 0,
        /// <summary>
        /// The value for the parameter will contain 1.0 in the Jacobian matrix.
        /// </summary>
        Frozen1 = 1,
        /// <summary>
        /// The value for the parameter will compute for the Jacobian matrix.
        /// </summary>
        Frozen2 = 2
    }

    /// <summary>
    /// This class provides the fitting of a curve by using the Gaussian least squares method.
    /// </summary>
    [Serializable]
    public class CurveFitting
    {
        /// <summary>
        /// Describes if the calculated curve have a precision error.
        /// </summary>
        private bool precisionError;

        /// <summary>
        /// The estimated error of the result.
        /// </summary>
        private double estimatedError;

        /// <summary>
        /// The needed iterations for fitting the curve.
        /// </summary>
        private int neededIterations;

        /// <summary>
        /// The model function which defines the curve to fit.
        /// </summary>
        private ModelFunction modelFunction;

        /// <summary>
        /// Initializes a new instance of the <see cref="CurveFitting"/> class.
        /// </summary>
        /// <param name="modelFunction">The model function which defines the curve to fit.</param>
        public CurveFitting(ModelFunction modelFunction)
        {
            this.modelFunction = modelFunction;
        }

        /// <summary>
        /// Describes if the fitting process have a precision error. If true the
        /// calculate curve parameters may not be correct.
        /// </summary>
        /// <value><c>True</c> if [precision error]; otherwise, <c>false</c>.</value>
        public bool PrecisionError
        {
            get { return this.precisionError; }
        }

        /// <summary>
        /// Gets the estimated error of the result.
        /// </summary>
        /// <value>The estimated error of the result.</value>
        public double EstimatedError
        {
            get { return this.estimatedError; }
        }

        /// <summary>
        /// Gets the needed iterations for fitting the curve.
        /// </summary>
        /// <value>The needed iterations for fitting the curve.</value>
        public int NeededIterations
        {
            get { return this.neededIterations; }
        }

        /// <summary>
        /// Gets or sets the model function which defines the curve to fit.
        /// </summary>
        /// <value>The model function which defines the curve to fit.</value>
        public ModelFunction ModelFunction
        {
            get { return this.modelFunction; }
            set { this.modelFunction = value; }
        }

        /// <summary>
        /// Fits the curve which is defined by the specified model function by using the Gaussian
        /// least squares method.
        /// </summary>
        /// <param name="x">The x values for fitting the curve.</param>
        /// <param name="y">The y values for fitting the curve.</param>
        /// <param name="parameters">The start values for the parameters to compute.</param>
        /// <param name="iterations">The maximum number of iterations to compute the parameters.</param>
        /// <returns>The computed parameters of the curve.</returns>
        public GeneralVector FitCurve(GeneralVector x, GeneralVector y, GeneralVector parameters, int iterations)
        {
            if (x == (GeneralVector) null)
            {
                throw new ArgumentNullException("x");
            }

            if (y == (GeneralVector) null)
            {
                throw new ArgumentNullException("y");
            }

            if (parameters == (GeneralVector) null)
            {
                throw new ArgumentNullException("parameters");
            }

            CurveFittingParameter[] parameterMask = new CurveFittingParameter[parameters.Count];

            for (int i = 0; i < parameterMask.Length; i++)
            {
                parameterMask[i] = 0;
            }

            return this.FitCurve(x, y, parameters, parameterMask, iterations, 1e-5);
        }

        /// <summary>
        /// Fits the curve which is defined by the specified model function by using the Gaussian
        /// least squares method.
        /// </summary>
        /// <param name="x">The x values for fitting the curve.</param>
        /// <param name="y">The y values for fitting the curve.</param>
        /// <param name="parameters">The start values for the parameters to compute.</param>
        /// <param name="iterations">The maximum number of iterations to compute the parameters.</param>
        /// <param name="precision">The precision of the parameters to compute.</param>
        /// <returns>The computed parameters of the curve.</returns>
        public GeneralVector FitCurve(GeneralVector x, GeneralVector y, GeneralVector parameters, int iterations,
                                      double precision)
        {
            CurveFittingParameter[] parameterMask = new CurveFittingParameter[parameters.Count];

            for (int i = 0; i < parameterMask.Length; i++)
            {
                parameterMask[i] = 0;
            }

            return this.FitCurve(x, y, parameters, parameterMask, iterations, precision);
        }

        /// <summary>
        /// Fits the curve which is defined by the specified model function by using the Gaussian
        /// least squares method.
        /// </summary>
        /// <param name="x">The x values for fitting the curve.</param>
        /// <param name="y">The y values for fitting the curve.</param>
        /// <param name="parameters">The start values for the parameters to compute.</param>
        /// <param name="parameterMask">The parameter mask, which defines the .</param>
        /// <param name="iterations">The maximum number of iterations to compute the parameters.</param>
        /// <returns>The computed parameters of the curve.</returns>
        public GeneralVector FitCurve(GeneralVector x, GeneralVector y, GeneralVector parameters,
                                      CurveFittingParameter[] parameterMask, int iterations)
        {
            return this.FitCurve(x, y, parameters, parameterMask, iterations, 1e-5);
        }

        /// <summary>
        /// Fits the curve which is defined by the specified model function by using the Gaussian 
        /// least squares method. 
        /// </summary>
        /// <param name="x">The x values for fitting the curve.</param>
        /// <param name="y">The y values for fitting the curve.</param>
        /// <param name="parameters">The start values for the parameters to compute.</param>
        /// <param name="parameterMask">The parameter mask, which defines the .</param>
        /// <param name="iterations">The maximum number of iterations to compute the parameters.</param>
        /// <param name="precision">The precision of the parameters to compute.</param>
        /// <returns>The computed parameters of the curve.</returns>
        public GeneralVector FitCurve(GeneralVector x, GeneralVector y, GeneralVector parameters,
                                      CurveFittingParameter[] parameterMask, int iterations, double precision)
        {
            if (x.Count != y.Count)
            {
                throw new ArgumentException("The number of x values has to be even with the number of y values.");
            }

            double eps = 0;
            int iterationCount = 0;
            int[] mask = new int[parameters.Count];
            GeneralVector b = new GeneralVector(parameters.Count);
            GeneralVector p = new GeneralVector(x.Count);
            double[,] pd = new double[x.Count,parameters.Count];
            double[,] pdsq = new double[parameters.Count,parameters.Count];

            for (int i = 0; i < parameters.Count; i++)
            {
                switch (parameterMask[i])
                {
                    case CurveFittingParameter.Free:
                        mask[i] = 0;
                        break;

                    case CurveFittingParameter.Frozen1:
                        mask[i] = 1;
                        break;

                    case CurveFittingParameter.Frozen2:
                        mask[i] = 2;
                        break;
                }
            }

            while (iterationCount < iterations)
            {
                double delta = 0;
                double tempuri = 0;
                GeneralVector r = this.modelFunction(x, parameters);

                for (int i = 0; i < parameters.Count; i++)
                {
                    if (mask[i] == 1)
                    {
                        for (int j = 0; j < x.Count; j++)
                        {
                            pd[j, i] = 1.0;
                        }
                    }
                    else
                    {
                        tempuri = parameters[i];
                        delta = (tempuri > 1.0) ? precision * tempuri : precision;
                        parameters[i] = tempuri + delta;
                        delta = parameters[i] - tempuri;
                        p = this.modelFunction(x, parameters);
                        parameters[i] = tempuri;

                        for (int j = 0; j < x.Count; j++)
                        {
                            pd[j, i] = (p[j] - r[j]) / delta;
                        }
                    }
                }

                for (int k = 0; k < parameters.Count; k++)
                {
                    for (int i = 0; i < parameters.Count; i++)
                    {
                        tempuri = 0.0;

                        for (int j = 0; j < x.Count; j++)
                        {
                            tempuri += pd[j, i] * pd[j, k];
                        }

                        pdsq[k, i] = tempuri;
                    }
                }

                if (!InvertMatrix(pdsq, parameters.Count))
                {
                    this.estimatedError = eps;
                    this.precisionError = true;
                    this.neededIterations = iterationCount;

                    return parameters;
                }

                for (int i = 0; i < x.Count; i++)
                {
                    p[i] = y[i] - r[i];
                    r[i] = p[i];
                }

                for (int i = 0; i < parameters.Count; i++)
                {
                    tempuri = 0.0;

                    for (int j = 0; j < x.Count; j++)
                    {
                        tempuri += pd[j, i] * p[j];
                    }

                    b[i] = tempuri;
                }

                for (int i = 0; i < parameters.Count; i++)
                {
                    tempuri = 0.0;

                    for (int j = 0; j < parameters.Count; j++)
                    {
                        tempuri += pdsq[i, j] * b[j];
                    }

                    p[i] = tempuri;
                }

                tempuri = 0.0;

                for (int i = 0; i < parameters.Count; i++)
                {
                    if (mask[i] != 0)
                    {
                        continue;
                    }

                    parameters[i] += p[i];
                    tempuri += p[i] * p[i];
                }

                eps = tempuri;

                if (tempuri < precision)
                {
                    break;
                }

                iterationCount++;
            }

            this.estimatedError = eps;
            this.precisionError = false;
            this.neededIterations = iterationCount;

            return parameters;
        }

        /// <summary>
        /// Inverts the specified matrix.
        /// </summary>
        /// <param name="a">The specified matrix.</param>
        /// <param name="n">The dimension of the matrix.</param>
        /// <returns>True, if the matrix is invertible otherwise, false.</returns>
        private static bool InvertMatrix(double[,] a, int n)
        {
            double[,] invertMatrix = new double[n,n];

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    invertMatrix[i, j] = 0.0;
                }

                invertMatrix[i, i] = 1.0;
            }

            for (int k = 0; k < n; k++)
            {
                double tempuri = a[k, k];

                if (tempuri == 0)
                {
                    return false;
                }

                for (int j = 0; j < n; j++)
                {
                    if (j > k)
                    {
                        a[k, j] /= tempuri;
                    }

                    invertMatrix[k, j] /= tempuri;
                }

                for (int i = 0; i < n; i++)
                {
                    if (i == k)
                    {
                        continue;
                    }

                    tempuri = a[i, k];

                    for (int j = 0; j < n; j++)
                    {
                        if (j > k)
                        {
                            a[i, j] -= a[k, j] * tempuri;
                        }

                        invertMatrix[i, j] -= invertMatrix[k, j] * tempuri;
                    }
                }
            }

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    a[i, j] = invertMatrix[i, j];
                }
            }

            return true;
        }
    }
}